{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting estimates out of DisMod-MR\n",
    "\n",
    "The goal of this document is to demonstrate how to export age-specific prevalence estimates from DisMod-MR in a comma-separated value (CSV) format, for use in subsequent analysis.\n",
    "\n",
    "It uses data from the replication dataset for regional estimates of HCV prevalence, as published in Mohd Hanafiah K, Groeger J, Flaxman AD, Wiersma ST. Global epidemiology of hepatitis C virus infection: New estimates of age-specific antibody to HCV seroprevalence. Hepatology. 2013 Apr;57(4):1333-42. doi: 10.1002/hep.26141. Epub 2013 Feb 4.  http://www.ncbi.nlm.nih.gov/pubmed/23172780\n",
    "\n",
    "The dataset is available from: http://ghdx.healthdata.org/sites/default/files/record-attached-files/IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
    "\n",
    "    wget http://ghdx.healthdata.org/sites/default/files/record-attached-files/IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
    "    unzip IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2019-07-23 13:51:38--  http://ghdx.healthdata.org/sites/default/files/record-attached-files/IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
      "Resolving ghdx.healthdata.org (ghdx.healthdata.org)... 2606:4700:10::6814:316, 2606:4700:10::6814:216, 104.20.3.22, ...\n",
      "Connecting to ghdx.healthdata.org (ghdx.healthdata.org)|2606:4700:10::6814:316|:80... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 74164 (72K) [application/zip]\n",
      "Saving to: 'IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP'\n",
      "\n",
      "100%[======================================>] 74,164      --.-K/s   in 0s      \n",
      "\n",
      "2019-07-23 13:51:38 (622 MB/s) - 'IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP' saved [74164/74164]\n",
      "\n",
      "Archive:  IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
      "  inflating: hcv_replication/2013_04_12_hcv_replication_code.ipynb  \n",
      "  inflating: hcv_replication/2013_04_12_hcv_replication_code.py  \n",
      "  inflating: hcv_replication/hierarchy.json  \n",
      "  inflating: hcv_replication/input_data.csv  \n",
      "  inflating: hcv_replication/nodes_to_fit.json  \n",
      "  inflating: hcv_replication/output_template.csv  \n",
      "  inflating: hcv_replication/parameters.json  \n",
      "  inflating: hcv_replication/README  \n",
      "  inflating: hcv_replication/README~  \n"
     ]
    }
   ],
   "source": [
    "!wget http://ghdx.healthdata.org/sites/default/files/record-attached-files/IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
    "!unzip IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This Python code will export predictions \n",
    "# for the following region/sex/year:\n",
    "predict_region = 'USA'\n",
    "predict_sex = 'male'\n",
    "predict_year = 2005"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import dismod code\n",
    "import dismod_mr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the model, and keep only data for the prediction region/sex/year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kept 20 rows of data\n"
     ]
    }
   ],
   "source": [
    "model_path = 'hcv_replication/'\n",
    "dm = dismod_mr.data.load(model_path)\n",
    "\n",
    "if predict_year == 2005:\n",
    "    dm.keep(areas=[predict_region], sexes=['total', predict_sex], start_year=1997)\n",
    "elif predict_year == 1990:\n",
    "    dm.keep(areas=[predict_region], sexes=['total', predict_sex], end_year=1997)\n",
    "else:\n",
    "    raise(Error, 'predict_year must equal 1990 or 2005')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using stored FE for beta_p_x_sex x_sex {'mu': 0.0, 'dist': 'Normal', 'sigma': 0.0001}\n",
      "finding initial values\n",
      ". . . \n",
      "finding MAP estimate\n",
      "\n",
      "finding step covariances estimate\n",
      "\n",
      "resetting initial values (1)\n",
      ". . . \n",
      "resetting initial values (2)\n",
      "\n",
      "mare: 0.36\n",
      "sampling from posterior\n",
      "\n",
      "CPU times: user 1min 22s, sys: 0 ns, total: 1min 22s\n",
      "Wall time: 1min 22s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<pymc.NormalApproximation.MAP at 0x7f6c0e29aba8>,\n",
       " <pymc.MCMC.MCMC at 0x7f6c15c746a0>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit model using the data subset (faster, but no borrowing strength)\n",
    "dm.vars += dismod_mr.model.process.age_specific_rate(dm, 'p', predict_region, predict_sex, predict_year)\n",
    "%time dismod_mr.fit.asr(dm, 'p', iter=2000, burn=1000, thin=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make posterior predictions\n",
    "pred = dismod_mr.model.covariates.predict_for(\n",
    "            dm, dm.parameters['p'],\n",
    "            predict_region, predict_sex, predict_year,\n",
    "            predict_region, predict_sex, predict_year, True, dm.vars['p'], 0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest way to get these predictions into a csv file is to use the Python Pandas package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This generates a csv with 1000 rows,\n",
    "# one for each draw from the posterior distribution\n",
    "\n",
    "# Each column corresponds to a one-year age group,\n",
    "# e.g. column 10 is prevalence at age 10\n",
    "\n",
    "pd.DataFrame(pred).to_csv(\n",
    "    model_path + '%s-%s-%s.csv'%(predict_region, predict_sex, predict_year))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 abie Domain Users 2.1M Jul 23 13:53 hcv_replication/USA-male-2005.csv\r\n"
     ]
    }
   ],
   "source": [
    "!ls -hal hcv_replication/$predict_region-*.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To aggregate this into pre-specified age categories, you need to specify the age weights and groups:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = [1, 8, 8, 9, 9, 10, 10, 10, 10, 10,\n",
    "           10, 10, 10, 10, 9, 9, 9, 9, 9, 9,\n",
    "           9, 9, 9, 9, 9, 9, 9, 9, 9, 9,\n",
    "           9, 9, 8, 8, 8, 8, 8, 8, 8, 8,\n",
    "           8, 7, 7, 7, 7, 7, 7, 7, 7, 7,\n",
    "           6, 6, 6, 6, 6, 6, 5, 5, 5, 5,\n",
    "           5, 5, 4, 4, 4, 4, 4, 4, 4, 3,\n",
    "           3, 3, 3, 3, 3, 3, 3, 3, 3, 3,\n",
    "           3, 1, 1, 1, 1, 1, 1, 1, 1, 1,\n",
    "           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAZxUlEQVR4nO3deZwlZX3v8c+XRZFFZRkIIjCoE71oAjoTBfUq0bhggrgLIRGXZJK4RF7GGDUxcb3RXJcrMaLjBUFjInjFMNcVLkIUF5BBGEBAUVFGiCAugLiBv/tHPV1zmOnpPtPT55yZns/79arXqXpq+9XT3efX9VTVU6kqJEkC2GbSAUiSNh8mBUlSz6QgSeqZFCRJPZOCJKm33aQD2BR77LFHLV68eNJhSNIWZdWqVT+oqkXTzduik8LixYu58MILJx2GJG1RknxnQ/NsPpIk9UwKkqSeSUGS1DMpSJJ6JgVJUs+kIEnqmRQkST2TgiSpZ1KQJPW26CeapdmsWrVqYvteunTpxPYtzZVnCpKknklBktQzKUiSeiYFSVLPpCBJ6pkUJEk9k4IkqedzChqLST4vIGl4nilIknomBUlSz6QgSeqZFCRJPZOCJKlnUpAk9UwKkqSeSUGS1BtZUkiyb5JzklyR5PIkL23lr03yvSQXt+FJA+u8KsnVSa5K8oRRxSZJmt4on2i+HfirqrooyS7AqiRntXnvqKq3Di6c5EDgKOCBwL2A/5fkN6vqjhHGKEkaMLIzhaq6vqouauO3AFcA+8ywypHAh6vqF1X1beBq4KGjik+StL6xXFNIshh4MHB+K3pxktVJTkqyayvbB7h2YLU1TJNEkixPcmGSC2+88cYRRi1JW5+RJ4UkOwMfBY6rqpuBE4D7AgcD1wNvm1p0mtVrvYKqFVW1rKqWLVq0aERRS9LWaaRJIcn2dAnhQ1V1OkBVfb+q7qiqXwPvY20T0Rpg34HV7w1cN8r4JEl3Nsq7jwKcCFxRVW8fKN97YLGnApe18ZXAUUnumuQAYAlwwajikyStb5R3Hz0C+GPg0iQXt7JXA0cnOZiuaega4M8AquryJKcBX6O7c+lF3nkkSeM1sqRQVecx/XWCT86wzpuAN40qJknSzHyiWZLUMylIknomBUlSz6QgSeqZFCRJPZOCJKlnUpAk9UwKkqSeSUGS1DMpSJJ6JgVJUs+kIEnqmRQkST2TgiSpZ1KQJPVMCpKknklBktQzKUiSeiYFSVLPpCBJ6pkUJEk9k4IkqWdSkCT1Zk0KSR40jkAkSZM3zJnCe5JckOSFSe458ogkSRMza1KoqkcCxwD7Ahcm+bckjxt5ZJKksRvqmkJVfQP4O+BvgEcDxye5MsnTRhmcJGm8hrmm8NtJ3gFcATwGOKKq/lsbf8eI45MkjdEwZwrvAi4CDqqqF1XVRQBVdR3d2cO0kuyb5JwkVyS5PMlLW/luSc5K8o32uWsrT5Ljk1ydZHWSh2z64UmSNsYwSeH0qvpgVf1sqmDqC76qPjjDercDf9XOKg4BXpTkQOCVwNlVtQQ4u00DHA4sacNy4ISNPRhJ0qYZJik8Z5qy5862UlVdP3BWcQtd89M+wJHAKW2xU4CntPEjgQ9U58vAPZPsPUR8kqR5st2GZiQ5GvhD4IAkKwdm7QLctDE7SbIYeDBwPrBXVV0PXeJIsmdbbB/g2oHV1rSy6zdmX5KkudtgUgC+SPeFvAfwtoHyW4DVw+4gyc7AR4HjqurmJBtcdJqymmZ7y+mal9hvv/2GDUOSNIQNJoWq+g7wHeDQuW48yfZ0CeFDVXV6K/5+kr3bWcLewA2tfA3dsxBT7g1cN01cK4AVAMuWLVsvaUiS5m6D1xSSnNc+b0ly88BwS5KbZ9twulOCE4ErqurtA7NWAse28WOBMwbKn9PuQjoE+MlUM5MkaTxmOlN4ZPvcZY7bfgTwx8ClSS5uZa8G3gycluQFwHeBZ7Z5nwSeBFwN3AY8b477lSTN0UwXmnebacWq+uEs889j+usEAI+dZvkCXjTTNiVJozXTheZVdBd6N3QB+D4jiUiSNDEzNR8dMM5AJEmTN1Pz0QOq6soNdTcx9WCaJGnhmKn56GV0zwO8bZp5RdchniRpAZmp+Wh5+/zd8YUjLRyrVq2ayH6XLl06kf1qYZjpTAGAJDsALwQeSXeG8HngPVX18xHHJkkas1mTAvABuq4t/rlNHw18kLXPF0iSFohhksL9q+qggelzklwyqoAkSZMzTNfZX23dTgCQ5GHAF0YXkiRpUma6JfVSumsI29P1SfTdNr0/8LXxhCdJGqeZmo/+YGxRSJI2C7N1nd1rL8PZYeQRSZImZtZrCkmenOQbwLeB/wSuAT414rgkSRMwzIXmNwCHAF9v/SE9Fi80S9KCNMwtqb+qqpuSbJNkm6o6J8lbRh6Z5t2knrCVtOUYJin8uL1n+fPAh5LcANw+2rAkSZMwTPPRkcDPgOOATwPfBI4YZVCSpMmY9Uyhqn6a5DeAhwI/BD5TVTeNPDJJ0tgNc/fRnwAXAE8DngF8OcnzRx2YJGn8hrmm8NfAg6fODpLsDnwROGmUgUmSxm+Yawpr6HpJnXILcO1owpEkTdJMfR+9rI1+Dzg/yRl0fR8dSdecJElaYGZqPtqlfX6zDVPOGF04kqRJmqnvo9cNTifZpSuuW0celSRpIoa5++hBSb4KXAZcnmRVkgeOPjRJ0rgNc6F5BfCyqtq/qvYH/gp432jDkiRNwjBJYaeqOmdqoqrOBXYaWUSSpIkZ5jmFbyV5DfDBNv1HdN1oS5IWmGHOFJ4PLAJOb8MewPNGGZQkaTJmTApJtgVeXVV/WVUPacNxVfWj2Tac5KQkNyS5bKDstUm+l+TiNjxpYN6rklyd5KokT9iko5IkzcmMSaGq7gCWznHbJwNPnKb8HVV1cBs+CZDkQOAo4IFtnXe3hCRJGqNhril8NclK4CPAT6cKq+r0mVaqqs8lWTxkHEcCH66qXwDfTnI1Xa+sXxpyfUnSPBjmmsJuwE3AY+jeo3AE8AebsM8XJ1ndmpd2bWX7cOf+lNa0svUkWZ7kwiQX3njjjZsQhiRpXUP1klpVP5in/Z1A987nap9vo7uQnWmWrek2UFUr6J6dYNmyZdMuI0mamw2eKSQ5IsmNwOoka5I8fFN3VlXfr6o7qurXdA/APbTNWgPsO7DovYHrNnV/kqSNM1Pz0ZuA/15V9wKeDvzjpu4syd4Dk0+l6zoDYCVwVJK7JjkAWII9sUrS2M3UfHR7VV0JUFXntw7xhpbk34HDgD2SrAH+ATgsycF0TUPXAH/Wtn95ktOArwG3Ay9qdz5JksZopqSw58A7Fdabrqq3z7Thqjp6muITZ1j+TXRnJ5KkCZkpKbyPte9UmG5akrTADP0+BUnSwjfMcwqSpK2ESUGS1JutQ7xtkjxrXMFIkiZrtg7xfg28eEyxSJImbJjmo7OSvDzJvkl2mxpGHpkkaeyG6fvo+e3zRQNlBdxn/sORJE3SrEmhqg4YRyCSpMmbNSkk2R74C+BRrehc4L1V9asRxiVJmoBhmo9OALYH3t2m/7iV/cmogpIkTcYwSeF3quqggenPJrlkVAFJkiZnmLuP7khy36mJJPcB7MFUkhagod68BpyT5Ft0b0jbH3jeSKOSJE3EMHcfnZ1kCXB/uqRwZVX9YuSRSZLGboNJIcljquqzSZ62zqz7JqGqTh9xbJKkMZvpTOHRwGeBI6aZV4BJQZIWmJnep/APSbYBPlVVp40xJknShNghniSpZ4d4kqSeHeJJknp2iCdJ6s3afJRkxyR/l2RFm16S5A9GH5okadyGuabwfuCXwMPb9BrgjSOLSJI0McMkhftW1T8BvwKoqp/RPdksSVpghkkKv0xyN7qLy7TO8ezmQpIWoGHuPnot8Glg3yQfAh6BHeJJ0oI0zN1HZyZZBRxC12z00qr6wcgjkySN3TB3H51dVTdV1Seq6uNV9YMkZw+x3klJbkhy2UDZbknOSvKN9rlrK0+S45NcnWR1kods2mFJkuZig0khyQ7tyeU9kuw68DTzYuBeQ2z7ZOCJ65S9Eji7qpYAZ7dpgMOBJW1YTve6T0nSmM3UfPRnwHF0CWAVa+84uhn4l9k2XFWfawlk0JHAYW38FOBc4G9a+QeqqoAvJ7lnkr2r6vqhjkKSNC9m6iX1ncA7k7ykqv55nva319QXfVVdn2TPVr4PcO3Acmta2XpJIclyurMJ9ttvv3kKS5IEw92S+l9JdgFoTzafPoI2/+mee6jpFqyqFVW1rKqWLVq0aJ7DkKSt2zBJ4TVVdUuSRwJPoGv2mWub//eT7A3QPm9o5WuAfQeWuzdw3Rz3IUmao2GSwh3t8/eBE6rqDOAuc9zfSuDYNn4scMZA+XPaXUiHAD/xeoIkjd8wD699L8l7gd8D3pLkrgx3K+u/011U3iPJGuAfgDcDpyV5AfBd4Jlt8U8CTwKuBm7Dh+MkaSKGSQrPoru19K1V9ePW7PPXs61UVUdvYNZjp1m2uPP7GiRJEzDrf/xVdRvwTeAJSV4M7FlVZ448MknS2A3TDPRS4EPAnm341yQvGXVgkqTxG6b56AXAw6rqpwBJ3gJ8CZivZxckSZuJYe4+CmvvQKKN+z4FSVqAhjlTeD9wfpKPtemnACeOLiRJ0qQM03X225OcCzyS7gzheVX11VEHJkkavw0mhSQ7AH8O3A+4FHh3Vd0+rsAkSeM30zWFU4BldAnhcOCtY4lIkjQxMzUfHVhVvwWQ5ETggvGEJEmalJnOFH41NWKzkSRtHWY6Uzgoyc1tPMDd2nToeqa4+8ijk7TRVq1aNZH9Ll26dCL71fya6SU7244zEEnS5A3z8JokaSthUpAk9UwKkqSeSUGS1DMpSJJ6JgVJUs+kIEnqmRQkST2TgiSpZ1KQJPVMCpKknklBktQzKUiSerO+o1nzb1JdG0vSbDxTkCT1TAqSpN5Emo+SXAPcAtwB3F5Vy5LsBpwKLAauAZ5VVT+aRHyStLWa5JnC71bVwVW1rE2/Eji7qpYAZ7dpSdIYbU7NR0cCp7TxU4CnTDAWSdoqTSopFHBmklVJlreyvarqeoD2ueeEYpOkrdakbkl9RFVdl2RP4KwkVw67YksiywH222+/UcUnSVuliZwpVNV17fMG4GPAQ4HvJ9kboH3esIF1V1TVsqpatmjRonGFLElbhbEnhSQ7Jdllahx4PHAZsBI4ti12LHDGuGOTpK3dJJqP9gI+lmRq//9WVZ9O8hXgtCQvAL4LPHMCsUnSVm3sSaGqvgUcNE35TcBjxx2PJGmtzemWVEnShJkUJEk9k4IkqWdSkCT1TAqSpJ5JQZLUMylIknomBUlSz6QgSeqZFCRJvUl1nT1xq1atmnQIkrTZ8UxBktQzKUiSeiYFSVJvq72mIGl+TfI63dKlSye274XGMwVJUs+kIEnqmRQkST2TgiSpZ1KQJPVMCpKknklBktQzKUiSeiYFSVLPpCBJ6tnNhaQt3qS62FiI3Wt4piBJ6pkUJEk9k4IkqbfZJYUkT0xyVZKrk7xy0vFI0tZks7rQnGRb4F+AxwFrgK8kWVlVX5tsZJK0voX4DonN7UzhocDVVfWtqvol8GHgyAnHJElbjc3qTAHYB7h2YHoN8LDBBZIsB5a3yVuTXDXHfe0B/GCO6y5E1sda1sVa1sWdLZT62H9DMza3pJBpyupOE1UrgBWbvKPkwqpatqnbWSisj7Wsi7WsizvbGupjc2s+WgPsOzB9b+C6CcUiSVudzS0pfAVYkuSAJHcBjgJWTjgmSdpqbFbNR1V1e5IXA58BtgVOqqrLR7S7TW6CWmCsj7Wsi7Wsiztb8PWRqpp9KUnSVmFzaz6SJE2QSUGS1FswSWG27jGS3DXJqW3++UkWt/Ldk5yT5NYk71pnnaVJLm3rHJ9kultmNzvzXRdJdkzyiSRXJrk8yZvHdzSbbhS/GwPrrkxy2WiPYP6M6O/kLklWJPl6+x15+niOZtOMqC6Obt8Zq5N8Oske4zmaeVRVW/xAd1H6m8B9gLsAlwAHrrPMC4H3tPGjgFPb+E7AI4E/B961zjoXAIfSPT/xKeDwSR/rJOoC2BH43TZ+F+DzW0JdjPJ3o81/GvBvwGWTPs5J1gXwOuCNbXwbYI9JH+sk6oLuxp0bpo4f+CfgtZM+1o0dFsqZwjDdYxwJnNLG/w/w2CSpqp9W1XnAzwcXTrI3cPeq+lJ1P+EPAE8Z6VHMj3mvi6q6rarOaeO/BC6ie4ZkSzDv9QGQZGfgZcAbRxf6vBtJXQDPB/4RoKp+XVVbwhO/o6iLtGGn1qpwd7bA56wWSlKYrnuMfTa0TFXdDvwE2H2Wba6ZZZubo1HURS/JPYEjgLM3OdLxGFV9vAF4G3Db/IQ5FvNeF+33AeANSS5K8pEke81fyCMz73VRVb8C/gK4lC4ZHAicOH8hj8dCSQqzdo8x5DKbsvzmYhR10a2UbAf8O3B8VX1rDrFNwrzXR5KDgftV1cc2JbAJGMXvxnZ0Z41fqKqHAF8C3jq38MZqFL8X29MlhQcD9wJWA6+aa4CTslCSwjDdY/TLtC+3ewA/nGWbg00kW0qXG6OoiykrgG9U1f+ahzjHZRT1cSiwNMk1wHnAbyY5d57iHaVR1MVNdGdLUwnyI8BD5iPYERtFXRwMUFXfbE3OpwEPn6+Ax2WhJIVhusdYCRzbxp8BfLb94KZVVdcDtyQ5pLUPPgc4Y/5Dn3fzXhcASd5I90dx3DzHO2qj+N04oaruVVWL6S44fr2qDpv3yOffKOqigP8LHNaKHgtsCe8/GcXfyfeAA5MsatOPA66Yx5jHY9JXuudrAJ4EfJ3ujoK/bWWvB57cxneg+y/marq7iu4zsO41dP8B3Er338GBrXwZcFnb5rtoT4Bv7sN81wXdf1FF9wt+cRv+ZNLHOcnfjYH5i9lC7j4aVV3QdcP8ObrmkrOB/SZ9nBOsiz9vfyer6ZLl7pM+zo0d7OZCktRbKM1HkqR5YFKQJPVMCpKknklBktQzKUiSeiYFzVmSpyapJA8Yw76OS7LjPG7vmk3pwTLJYUk+3safPF0vm3PY5rlJ1nspfCu/KsklSb6Q5P6buq9Z4nhtkpePch/afJkUtCmOpnui96gx7Os4ut5aJyLJthuaV1Urq2rU3YkfU1UH0XXQ9j/XnTlTfNLGMCloTlovoY8AXsBAUkiyTZJ3t/cufDzJJ5M8o81bmuQ/k6xK8pnWE+26290p3bsbLklyWZJnJ/lLur5kzklyTlvuhCQXtv28bmD9a5K8rnXOdunUWUzrA//MJF9N8l4G+rVJ8h8tpsuTLB8ovzXJ65OcDxyarv/9K5OcR9dt9tRyz53qVz/JxQPDz5I8uh3TSUm+0vZ/ZFv2bkk+nK7v/VOBuw1R9Z8D7jdwrH/f4nlmkvum68N/VZLPJ3lAknu05bZp6+yY5Nok2yf50xbTJUk+Ot2Z2HTbbOUnp3vHyBeTfGvqZ9zmvaLV/SVp797Y0Ha0GZr003MOW+YA/BFwYhv/IvCQNv4M4JN0/3D8BvCjVrZ9W25RW+7ZwEnTbPfpwPsGpu/RPq9hoJ9+YLf2uS1wLvDbA8u9pI2/EPjfbfx44O/b+O/TPaG9xzrbuhvdE+y7t+kCntXGd6DrMXMJXUI5Dfh4m/dc1n/HwBF0753YHvgfwB+18nvSPUW7E13X2ye18t8GbgeWTVMn506VA3/N2n79rwFeMbDc2cCSNv4wum4ZoOueZep9GM8eqJPdB9Z940C9vRZ4+SzbPJnuad9t6J56v7qVH95+zjuuU7fTbsdh8xu2Q5qbo4GpjvE+3KYvousL6CNV9Wvgv6b+swfuDzwIOCvdC+y2Ba6fZruXAm9N8ha6L93Pb2D/z2r/1W8H7E33xbS6zTu9fa5i7X/0j5oar6pPJPnRwLb+MslT2/i+dF/8NwF3AB9t5Q8Avl1V3wBI8q/AcqaRZAldE89jqupXSR4PPHmgnX4HYL8W0/EtptVJVk+3veZDSX5GS3oD5ae2fe5M1/naR7L2BYF3HVjm2cA5dGd1727lD0rXp9U9gZ2Bz6xzHDNtE+A/2s/5a1nbXfbvAe+vqtvacf1wiO1oM2JS0EZLsjvwGLovlaL7gq8kr2D67oZp5ZdX1aHrbGtfuj5ioHvL1XuSLKXrl+Yfk5xZVa9fZ50DgJcDv1NVP0pyMt0X7ZRftM87uPPv+Hp9uiQ5jO6L7NCqui1db6dT2/p5Vd0x0/rTbG8nurOIP62qqV43Azy9qq5aZ9mhttkcU1UXTlP+0/a5DfDjqjp4mmVW0tXlbsBS4LOt/GTgKVV1SZLnsrZTuykzbRPW1jOs/bmH9Y9ptu1oM+I1Bc3FM4APVNX+VbW4qvYFvk13lnAe8PR01xb2Yu0XzVXAoiSHQtf3fJIHVtW1VXVwG96T5F7AbVX1r3T98k91w3wLsEsbvzvdl+FP2j4OHyLmzwHHtH0fDuzayu8B/KglhAcAh2xg/SuBA5Lct00fvYHl3k/3n/LgGc5ngJekZYEkD54mpgfRNSHNSVXdDHw7yTPb9pLkoDbvVroO3d5Jd/Y1leh2Aa5P9x6AYzZmmzM4E3j+1PWJJLvNcTuaEJOC5uJo1vafP+WjwB+2zzV0bfPvBc4HflLdKw+fAbwlySV0Pa1O19f8bwEXJLkY+FvWvu5yBfCpJOdU1SXAV4HLgZOALwwR8+uARyW5CHg88N1W/mlgu9Z08wbgy9OtXFU/p2su+kS7sPuddZdJsn87xudn7cXmZW272wOrk1zWpgFOAHZu+34F3Rf3pjgGeEGr38u58+slT6W7DnTqQNlr6H4+Z9ElvY3d5nqq6tN0ZyYXtp/hVJPZRm1Hk2MvqZp3SXauqltbM9MFwCOq6r8mHZek2XlNQaPw8XTv7r0L8AYTgrTl8ExBktTzmoIkqWdSkCT1TAqSpJ5JQZLUMylIknr/HwqQIaxK3bFSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 1000 samples from the posterior distribution for age-standardized prevalence\n",
    "import numpy as np, matplotlib.pyplot as plt\n",
    "\n",
    "age_std = np.dot(pred, weights) / np.sum(weights)\n",
    "plt.hist(age_std, color='#cccccc', density=True)\n",
    "plt.xlabel('Age-standardized Prevalence')\n",
    "plt.ylabel('Posterior Probability');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can extract an age-standardized point and interval estimate from the 1000 draws from the posterior distribution stored in age_std as follows:  (to do this for crude prevalence, use the population weights to average, instead of standard weights.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "age_std prev mean: 0.013684724543190618\n",
      "age_std prev 95% UI: [0.01056258 0.01612552]\n"
     ]
    }
   ],
   "source": [
    "import pymc as mc\n",
    "\n",
    "print('age_std prev mean:', age_std.mean())\n",
    "print('age_std prev 95% UI:', mc.utils.hpd(age_std, .05))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For groups, just do the same thing group by group:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_cutpoints = [0, 1,  5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   a0  a1   mu  std\n",
      "0   0   1  0.0  0.0\n",
      "1   1   5  0.0  0.0\n",
      "2   5  10  0.0  0.0\n",
      "3  10  15  0.0  0.0\n",
      "4  15  20  0.0  0.0\n"
     ]
    }
   ],
   "source": [
    "results = []\n",
    "for a0, a1 in zip(group_cutpoints[:-1], group_cutpoints[1:]):\n",
    "    age_grp = np.dot(pred[:, a0:a1], weights[a0:a1]) / np.sum(weights[a0:a1])\n",
    "    results.append(dict(a0=a0,a1=a1,mu=age_grp.mean(),std=age_grp.std()))\n",
    "    \n",
    "results = pd.DataFrame(results)\n",
    "print(np.round(results.head(), 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAASqklEQVR4nO3df4xlZ13H8ffHXbujGFtZFqO7rbukC7ooRXcsqEg2VMzWHyyJbdhqYv9oshpp/B0tMTSl4Q9qjBVDgza0WjaGFssPJ7jaKGU1GqjdFSwstTKUaseiLGwtFrmUha9/3DN6uXtn58zMvTN37n2/kpu555znzDxnzsz9nOc55zwnVYUkabp93UZXQJK08QwDSZJhIEkyDCRJGAaSJGDrRleg33Oe85zavXv3RldDkjaVkydPfraqdqx2/bELg927d3PixImNroYkbSpJ/nUt69tNJEkyDCRJhoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANp0zhw4AAHDhzY6GpoQhkGkiTDQGrLI3NNMsNAktQuDJIcTPJIkvkkNwxYvi3JPc3yB5LsbubvTvLFJB9pXn8w3OpLkoZh2SGsk2wBbgNeCSwADyaZq6qP9xS7Dniyqi5Nchi4BXhNs+yTVfXiIddbkjREbVoGlwPzVfVoVT0D3A0c6itzCLireX8vcEWSDK+akqRRahMGO4HHe6YXmnkDy1TVWeApYHuzbE+SDyf5myQ/vMb6SpJGoM2TzgYd4VfLMp8GLqmqzyXZD7w3yQur6vNfs3JyBDgCcMkll7SokiRpmNq0DBaAi3umdwFPLFUmyVbgQuBMVX2pqj4HUFUngU8Cz+//AVV1e1XNVtXsjh2rfoSnJGmV2oTBg8DeJHuSXAAcBub6yswB1zbvrwLur6pKsqM5AU2S5wF7gUeHU3VJ0rAs201UVWeTXA/cB2wB7qyqU0luBk5U1RxwB3A0yTxwhm5gALwcuDnJWeArwM9X1ZlRbIgkafXanDOgqo4Bx/rm3djzvgNcPWC9dwHvWmMdJUkj5h3I0jpxOAuNM8NAkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMpKnhTW86H8NAkmQYSJIMA0kShoEkCcNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEoaBJAnDQJKEYSBJwjCQJGEYSJIwDCRJGAaSJAwDSRKGgaQWfH7y5DMMJEntwiDJwSSPJJlPcsOA5duS3NMsfyDJ7r7llyR5OsmvD6fakqRhWjYMkmwBbgOuBPYB1yTZ11fsOuDJqroUuBW4pW/5rcBfrL26kqRRaNMyuByYr6pHq+oZ4G7gUF+ZQ8Bdzft7gSuSBCDJq4FHgVPDqbIkadi2tiizE3i8Z3oBeMlSZarqbJKngO1Jvgj8JvBKYMkuoiRHgCMAl1xySevKS5tBp9NhZmaG48ePD5wvjYM2LYMMmFcty7wBuLWqnj7fD6iq26tqtqpmd+zY0aJK0uYxMzNDknNeBoHGSZuWwQJwcc/0LuCJJcosJNkKXAicoduCuCrJbwMXAV9N0qmqt6y55tIYGnT5ZX+L4HzlB5W1ZaH10CYMHgT2JtkD/DtwGPjpvjJzwLXAB4GrgPurqoAfXiyQ5CbgaYNAWpnFlkW/7r+YNBzLhkFzDuB64D5gC3BnVZ1KcjNwoqrmgDuAo0nm6bYIDo+y0tJqLB6Fn+9Ifa1W+r37y4+iZSG10aZlQFUdA471zbux530HuHqZ73HTKuonbTi7aTQNWoWBNM3W2k3T6XQGlh0UJmttWUirZRhIPUbRTbP4gd/fTWWrQuPEsYmkMbfYsuh/dTqdja6aJogtA6nHOHbT2LLQerBlIEkyDKTl2E2jaWAYSMvo7abpPWFsN40miWEgSTIMJI2ej80cf4aBJMkwkCQZBpIkvOlMmngOtKc2DANpwvk8BLVhGEgTxOchaLU8ZyBJsmUgTZJxHGhPm4MtA0mSYSBNOgfaUxuGgTThHGhPbRgGkiTDQJJkGEiSMAwkSRgGkjYBn4cweoaBJMkwkCQ5HIWk83D46+lhGEhaksNfTw/DQNL/WcmQ1g5/PVk8ZyBJatcySHIQeDOwBXhbVb2pb/k24O3AfuBzwGuq6rEklwO3LxYDbqqq9wyr8pKGayVH97YEJsuyLYMkW4DbgCuBfcA1Sfb1FbsOeLKqLgVuBW5p5n8MmK2qFwMHgT9MYteUJI2ZNh/MlwPzVfUoQJK7gUPAx3vKHAJuat7fC7wlSarqf3rKzACeddLU2oxH0ovDXw+a79VEk6XNOYOdwOM90wvNvIFlquos8BSwHSDJS5KcAj4K/Hyz/GskOZLkRJITp0+fXvlWSBoJh7+eHm3C4Nzrys49wl+yTFU9UFUvBL4feF2Sc/6Kqur2qpqtqtkdO3a0qJIkaZjahMECcHHP9C7giaXKNOcELgTO9BaoqoeBLwDfvdrKSpJGo00YPAjsTbInyQXAYWCur8wccG3z/irg/qqqZp2tAEm+A3gB8NhQai5JGpplTyBX1dkk1wP30b209M6qOpXkZuBEVc0BdwBHk8zTbREcblZ/GXBDki8DXwV+oao+O4oNkSStXqvLPKvqGHCsb96NPe87wNUD1jsKHF1jHSVJI+Y1/5p4DrYmLc8w0MRzsDVpeYaBJsqgp2E52Jq0PMNA0sjYRbd5GAaaKCs9srclMFp20W0ehoGkobCLbnPzeQaaeIuDrfW/Op3ORldNGhu2DDTxegdbg/8/ArXPerjsotvcbBlIkgwDSaNjF93mYRhIGhmfh7B5GAaSJE8gS5uFJ1w1SrYMJGmN+rvBNiNbBlJLHplrktkykCQZBpIkw0CShOcMJI0xh8BeP4aBpLHlENjrxzCQNBYcAntjGQbSlPDDUudjGEgaCw6BvbG8mkiSZBhIGl8Ogb1+DANJY8shsNePYSBJMgwkSYaBJAnDQJKEYSBJomUYJDmY5JEk80luGLB8W5J7muUPJNndzH9lkpNJPtp8fcVwqy9pPRw/ftybvCbcsmGQZAtwG3AlsA+4Jsm+vmLXAU9W1aXArcAtzfzPAj9ZVd8DXAscHVbFJUnD06ZlcDkwX1WPVtUzwN3Aob4yh4C7mvf3AlckSVV9uKqeaOafAmaSbBtGxSVJw9MmDHYCj/dMLzTzBpapqrPAU8D2vjI/BXy4qr7U/wOSHElyIsmJ06dPt627JGlI2oTBuYOJQ/9g4uctk+SFdLuOfm7QD6iq26tqtqpmd+zY0aJKkqRhahMGC8DFPdO7gCeWKpNkK3AhcKaZ3gW8B/jZqvrkWissSRq+NmHwILA3yZ4kFwCHgbm+MnN0TxADXAXcX1WV5CLgz4HXVdXfD6vSkqThWjYMmnMA1wP3AQ8D76yqU0luTvKqptgdwPYk88CvAouXn14PXAq8PslHmtdzh74VkqQ1afVwm6o6Bhzrm3djz/sOcPWA9d4IvHGNdZQkjZh3IEuSDANJkmEgSaLlOQNpHCw+6coxcjYf99n4s2UgSTIMJEmGgSQJzxloithvLS3NloEkyTCQJBkGkiQMA0kShoEkCcNAkoSXlkrSqnU6HWZmZs65bHlx/mZiGEjSKs3MzJCc+wj4qv7HxI8/w0DSRFgcyHA9ne9Gxo2oz1p4zkCSZMtA0vhrM5TIuA03st71GdRdtRK2DCRJhoEkrVan06Gqznl1Op2NrtqKGQaStEqLl48eOHDga04Yb7bLSsEwkCThCWRtApN0Y480rgwDjb1JurFHGleGgcbKoBt1VnJjz7hdXihtFp4zkCTZMtB4WemRvS0BaThsGUiSDAONv0m6sUcaV4aBxt4k3dgjjatWYZDkYJJHkswnuWHA8m1J7mmWP5BkdzN/e5IPJHk6yVuGW3VJ0rAsGwZJtgC3AVcC+4BrkuzrK3Yd8GRVXQrcCtzSzO8Arwd+fWg1liQNXZuWweXAfFU9WlXPAHcDh/rKHALuat7fC1yRJFX1har6O7qhIEkaU23CYCfweM/0QjNvYJmqOgs8BWxvW4kkR5KcSHLi9OnTbVeTJA1JmzAY9MSE/nEA2pRZUlXdXlWzVTW7Y8eOtqtJkoakTRgsABf3TO8CnliqTJKtwIXAmWFUUJOj/2ogSeOjTRg8COxNsifJBcBhYK6vzBxwbfP+KuD+chQxSdo0lh2OoqrOJrkeuA/YAtxZVaeS3AycqKo54A7gaJJ5ui2Cw4vrJ3kM+GbggiSvBn60qj4+/E2RJK1Wq7GJquoYcKxv3o097zvA1Uusu3sN9ZMkrQPvQJYkGQaSJMNAkoRhIEnCMJAkYRhIkjAMJEkYBpIkDANJEi3vQJbGwfHjxze6CtLEsmUgSTIMJEmGgSQJw2BdrfXhLj4cRtKoGAaSJK8m0uh1Oh1mZmbOuRpocb6kjWcYaORmZmZIcs58n4wqjQ/DQEM16JzG+e4P6C/vvQTSxvCcgSTJloGGa6VH9rYEpPFgGEjSGk3CQY3dRGpttfc5dDodquqcV6fTGX4lJa2KYbAC3vS1OouXj/b//rysVBofhoEkyTCYJrZsJC3FE8halncQS5PPMNCyvINYmnyGwTpY65H1eh6ZewexNJ2mMgwWP8DafnCt9cN4rUfWHplLGrWpDIOVWsmH8VqPrAdZ65H5SsLMO4il6TRVYdDmQ3EUH+YbzZaFpOW0CoMkB4E3A1uAt1XVm/qWbwPeDuwHPge8pqoea5a9DrgO+Arwi1V139Bq32epbpvF+evxobjeR9b964+iz3/xDuJB872aSJoMy4ZBki3AbcArgQXgwSRzVfXxnmLXAU9W1aVJDgO3AK9Jsg84DLwQ+Hbgr5M8v6q+MuwNgfMfAR84cKDVh6LdJOfqbzUtbrNBIE2ONi2Dy4H5qnoUIMndwCGgNwwOATc17+8F3pLup/Ih4O6q+hLwqSTzzff74HCqvzms9ch6JesbZpJWo00Y7AQe75leAF6yVJmqOpvkKWB7M/9Dfevu7P8BSY4AR5rJp5M80qr2ffbv379/0PzZ2VlOnjx5cnZ29rzLm7qcs/xFL3rRZfv37z/nd3XZZZedfeihh/5pNXUdoecAn+2dsdzvZaU/YNDvaD3XX8Y52z9lpnn7p3nbAV6wlpXbhMGg/9z+w9SlyrRZl6q6Hbi9RV2WleREVc0O43ttRm6/2z+t2z/N2w7d7V/L+m3GJloALu6Z3gU8sVSZJFuBC4EzLdeVJG2wNmHwILA3yZ4kF9A9ITzXV2YOuLZ5fxVwf3U7ueeAw0m2JdkD7AX+YThVlyQNy7LdRM05gOuB++heWnpnVZ1KcjNwoqrmgDuAo80J4jN0A4Om3Dvpnmw+C7x2VFcS9RhKd9Mm5vZPt2ne/mnedljj9scbjyRJPs9AkmQYSJImLAySHEzySJL5JDdsdH1GKcnFST6Q5OEkp5L8UjP/2Un+Ksknmq/fstF1HaUkW5J8OMn7muk9SR5otv+e5qKHiZTkoiT3Jvnn5u/gB6Zp/yf5leZv/2NJ3pFkZpL3f5I7k3wmycd65g3c3+n6/eaz8KEk37fc95+YMOgZNuNKYB9wTTMcxqQ6C/xaVX0X8FLgtc323gC8v6r2Au9vpifZLwEP90zfAtzabP+TdIdKmVRvBv6yqr4TuIzu72Eq9n+SncAvArNV9d10L25ZHApnUvf/HwMH++Yttb+vpHv15l66N/S+dblvPjFhQM+wGVX1DLA4bMZEqqpPV9U/Nu//m+4HwU6623xXU+wu4NUbU8PRS7IL+HHgbc10gFfQHRIFJnj7k3wz8HK6V/JRVc9U1X8xRfuf7tWQ39Dc2/SNwKeZ4P1fVX9L92rNXkvt70PA26vrQ8BFSb7tfN9/ksJg0LAZ5wx9MYmS7Aa+F3gA+Naq+jR0AwN47sbVbOR+D/gN4KvN9Hbgv6rqbDM9yX8DzwNOA3/UdJO9LcmzmJL9X1X/DvwO8G90Q+Ap4CTTs/8XLbW/V/x5OElh0Groi0mT5JuAdwG/XFWf3+j6rJckPwF8pqp6x1aapr+BrcD3AW+tqu8FvsCEdgkN0vSNHwL20B0R+Vl0u0b6Ter+X86K/xcmKQymbuiLJF9PNwj+pKre3cz+z8XmYPP1MxtVvxH7IeBVSR6j2yX4CrothYuabgOY7L+BBWChqh5opu+lGw7Tsv9/BPhUVZ2uqi8D7wZ+kOnZ/4uW2t8r/jycpDBoM2zGxGj6x+8AHq6q3+1Z1Ds0yLXAn6133dZDVb2uqnZV1W66+/r+qvoZ4AN0h0SByd7+/wAeT7I4UuUVdO/0n4r9T7d76KVJvrH5X1jc/qnY/z2W2t9zwM82VxW9FHhqsTtpSVU1MS/gx4B/AT4J/NZG12fE2/oyus2+h4CPNK8fo9tv/n7gE83XZ290Xdfhd3EAeF/z/nl0x7+aB/4U2LbR9Rvhdr8YONH8DbwX+JZp2v/AG4B/Bj4GHAW2TfL+B95B9/zIl+ke+V+31P6m2010W/NZ+FG6V12d9/s7HIUkaaK6iSRJq2QYSJIMA0mSYSBJwjCQJGEYSJIwDCRJwP8CQbz4Tu36T4IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(.5*(results.a0+results.a1), results.mu,\n",
    "             xerr=.5*(results.a1-results.a0),\n",
    "             yerr=1.96*results['std'],\n",
    "             fmt='ks', capsize=0, mec='w')\n",
    "plt.axis(ymin=0, xmax=100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm IHME_GBD_HEP_C_RESEARCH_ARCHIVE_Y2013M04D12.ZIP\n",
    "!rm -r hcv_replication/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tue Jul 23 13:53:06 PDT 2019\r\n"
     ]
    }
   ],
   "source": [
    "!date"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dismod_mr",
   "language": "python",
   "name": "dismod_mr"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
